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Abstract 

The direct simulation Monte Carlo (DSMC) method is the established technique for 
the simulation of rarefied gas flows. In some flows of engineering interest, such as occur 
for aero-braking spacecraft in the upper atmosphere, DSMC can become prohibitively 
expensive in CPU time because some regions of the flow, particularly on the windward 
side of blunt bodies, become collision dominated. As an alternative to using a hybrid 
DSMC and continuum gas solver (Euler or Navier-Stokes solver) this work is aimed 
at making the particle simulation method efficient in the high density regions of the 
flow. A high density, infinite collision rate limit of DSMC, the Equilibrium Particle 
Simulation method (EPSM) was proposed some 15 years ago. EPSM is developed here 
for the flow of a gas consisting of many different species of molecules and is shown to 
be computationally efficient (compared to DSMC) for high collision rate flows. It thus 
offers great potential as part of a hybrid DSMC/EPSM code which could handle flows 
in the transition regime between rarefied gas flows and fully continuum flows. As a 
first step towards this goal a pure EPSM code is described. The next step of combining 
DSMC and EPSM is not attempted here but should be straightforward. EPSM and 
DSMC are applied to Taylor- Couette flow with Kn = 0.02 and 0.0133 and S w = 3). 
Toroidal vortices develop for both methods but some differences are found, as might 
be expected for the given flow conditions. EPSM appears to be less sensitive to the 
sequence of random numbers used in the simulation than is DSMC and may also be 
more dissipative. The question of the origin and the magnitude of the dissipation in 
EPSM is addressed. It is suggested that this analysis is also relevant to DSMC when 
the usual accuracy requirements on the cell size and decoupling time step are relaxed 
in the interests of computational efficiency. 


*This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-19480 while the author was in residence at the Institute for Computer Applications for 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


The Direct Simulation Monte Carlo (DSMC) method [1] [2] is widely used for high speed 
rarefied gas flows for which the continuum equations of fluid mechanics are expected to be 
invalid. In this flow regime, the continuum gas properties such as density, flow velocity and 
temperature must be supplemented with information about the distribution of molecular 
energies amongst the molecules of the gas. In the case that the only relevant energies are 
those of translation and rotation of the molecules, the distribution function /(v, e rot , r, t) is 
defined such that, at a given point r in the gas, at a given instant t, the expected number 
density (number/ volume) of molecules which have a velocity in the range v to v + dv and 
also have a value of rotational energy in the range £ rot to e TOt + de ro t is given by nfdvde TOt , 
where n is the number density of all molecules at the point. 

In DSMC the molecules of a gas are represented by a very much smaller number of 
‘simulator particles’. The particles are first moved in collisionless trajectories for a small 
time interval At; then the particle velocities are subject to alteration during the collision 
phase of the calculation. The flow field is divided into small cells and representative collisions 
are calculated amongst the particles in each cell while conserving energy and momentum in 
each collision. The number of collisions calculated in each cell in each time step At reflects 
the appropriate local collision rate in the real gas. If the collision rate is large enough the 
distribution of the finite sample of particle velocities in a cell will conform (in a statistical 
sense) to the local Maxwell- Boltzmann equilibrium distribution / 0 (see equation (1)). If 
the collision rate is low there will be little re-distribution of molecular velocities and / will 
remain close to the value established in the cell by the free flight part of the simulation, and 
can thus be far removed from f 0 . 

In many blunt body atmospheric re-entry flows the gas near the windward surfaces of 
a spacecraft can be at a high density and the continuum fluid dynamics equations should 
apply, while the flow near the leeward surfaces is at a low density where DSMC should be 
used. It is appealing to use a particle method throughout these simulations rather than a 
hybrid DSMC/continuum code in which, to mention some apparent difficulties, the inter- 
face between the continuum and particle solvers could be bothersome [3] and the continuum 
solver might be adversely affected by the statistical noise generated by the DSMC calcula- 
tions. The flow in the high density regions is collision dominated and, when using a particle 
method, prohibitively expensive amounts of CPU time can be spent calculating collisions. A 
suggestion by Bird [4] to limit the number of collisions in any cell in any time step to (say) 
10 per particle has sometimes been used. The justification for this is that after 10 collisions 
per particle have been calculated the state in the cell will, in effect, be at equilibrium and 
more collisions cannot improve the accuracy of the simulation. However, even this collision 
limiting method can be prohibitively expensive. 

The purpose of this work is to describe a high density version of DSMC which will be 
suitable for later use as part of a hybrid particle simulation code. The high density direct 
simulation method is called the Equilibrium Particle Simulation Method (EPSM) and was 
proposed by Pullin [5] as an infinite collision rate limit of DSMC for a single species of particle 
only but it contains all the necessary features of a method which can be easily combined 
with DSMC. In a hybrid DSMC/EPSM code the EPSM routines can be used to establish 
a new distribution of energies amongst the particles in those cells for which the theoretical 
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collision rate would require more than (say) 10 collisions per particle and this can be done 
without calculating any collisions at all in those cells. The two developments of EPSM which 
must be made before it can be suitable for a hybrid DSMC/EPSM code are 

1. rotation energy must be handled in a manner compatible with existing DSMC codes 
and 

2. multiple species of particles must be allowed for. 

These two issues are addressed in sections (2) and (3). 

A new code was constructed by replacing the collision routines in a freely available DSMC 
code with EPSM routines. The original code, called DSMC2A, was written by G. A. Bird 
[2] and can be obtained from the computer disk which accompanies ref. [2]. DSMC2A deals 
with a 2D axisymmetric flow in a rectangular region and has various options for setting 
boundary conditions and placing a surface within the flow; it allows for different species of 
molecules, any of which may possess energy of molecular rotation and vibration as well as 
translational energy. Some small coding errors in the original DSMC2A code were corrected. 
These are listed in Appendix B. 

The new EPSM code is tested for Taylor-Couette flow and the results compared with the 
DSMC results obtained by repeating the DSMC calculations made at low Knudsen number 
by Stefanov and Cercignani [6] with the same grid size. It has been found that similar, but 
not identical, vortices will form using EPSM for which there is no attempt to model the 
collision rate (and hence the dissipation) correctly. The origin of the dissipation in EPSM is 
discussed in section 7; that discussion has some implications for DSMC when it is applied 
in the near continuum regime. 


2 Vibration energy in EPSM 


The Equilibrium Particle Simulation Method (EPSM) is similar to DSMC except that there 
is no calculation of molecular collisions. Instead, where DSMC would calculate collisions, it 
is assumed in EPSM that the collision rate is always large enough to bring about within a 
cell a distribution of particle velocities which represents (within the limitations set by the 
finite sample size) an approximation to an equilibrium distribution. For a single species of 
particle, with a particle mass m s , for which rotation of the molecules (with two degrees of 
freedom) is active, the equilibrium distribution function is 


/o(v, Srot) 


f— ) 

\2irhTj 


exp 


— m s (v- v) 2 \ 
2 k b T ) 



( 1 ) 


where kb is Boltzmann’s constant. The distribution of a single component of velocity, say 
v x , can be obtained from (1) by performing the integration / 0 °° fodv y dv z de r ot- The 

result is 



( 2 ) 


Thus a single component of molecular velocity is normally distributed in an equilibrium gas. 
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EPSM bypasses the calculation of a large number of collisions amongst the finite sample 
of particles in each cell by selecting the new velocity components from a normal distribution 
in such a way that the mean and variance of the finite sample conform to values determined 
by the requirement that the total energy and momentum be conserved in each cell. This 
problem is different from the simpler problem of generating a finite sample of velocities 
selected from a parent distribution which has a known mean and variance. That is, because 
of the finite sample of simulator particles which are used in any simulation (or any ensemble 
of simulations) the values of v and T in (1) are not known exactly in EPSM or DSMC; 
what is known is the estimate of these quantities which are derived from the finite sample. 
Thus the estimate of the macroscopic fluid velocity, the mass average particle velocity v est , 
is given by 

N 

Vest = Y mkVk/Mtot (3) 

k = 1 

where N is the number of particles in the sample, is the mass of the k-th particle in the 
sample and 

N 

Mtot = Y mk ( 4 ) 

fc=l 

is the total mass of the N particles. The estimate of the temperature, T est , is derived from 
the total energy of the finite sample of N particles. It is related to the average thermal 
energy per single degree of freedom for a single particle, e' d o j (see equation (11)), by 


& d.o.f — cykbT es t- 


( 5 ) 


In the original application of EPSM [5] a one dimensional flow of a single species of particle 
was treated. The mathematical derivation of a procedure for selecting velocity components 
which are normally distributed and are constrained to have a specified mean and variance 
was developed in [7]. The procedure is given in the form of an algorithm in [5]; the algorithm 
is also given here in Appendix A but in greater detail, particularly for the special cases when 
the number of velocities to be set is low. 

In the original application of EPSM [5] energy of molecular structure, such as rotation, 
was not assigned to individual particles in the simulation but rather the total energy of 
molecular structure was calculated by assuming each degree of freedom contained exactly the 
same thermal energy as was contained in a single translational degree of freedom^. In current 
DSMC codes energy of molecular rotation is assigned to each particle, and for compatibility 
with these codes the EPSM code developed here allows for molecular rotation energy as 
described below. The code can be easily generalized to any number of fully excited (classical) 
degrees of freedom but the issue of treating vibration energy which is not fully excited is not 
dealt with here. 

It is convenient to represent molecular rotation in two degrees of freedom by a pseudo- 
velocity having two components (/j, i,^ 2 ), where 


Srot — T ^ • 


(6) 


tin fact a similar argument was applied to the energy stored in the two redundant velocity components 
which then did not have to be stored for the one dimensional flow case. 
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Then, in an equilibrium state, the distribution function is given by 


fo(v x ,v y ,v z ,fii,n 2 ) 



ft 2 - / 2 , t 2 , 2 

X +Uy + t> 2 + /Xj 


2 k b T 



(V 


where v' x = v x — v x , v' y = v y — v yi v' z = v z — v z are the three components of thermal velocity. 
If an approximation to the equilibrium distribution (7) can be established in a cell the two 
components of pseudo-velocity can be converted to a single value of rotational energy e TOt 
which is stored for each particle. It follows from (6) and (7) that the distribution of e TOt so 
derived will conform to (1). 


3 Multiple species 


The EPSM code used here was derived by replacing the collision routines in DSMC2A 
(subroutines COLLMR, SELECT2A, ELASTIC, INELR and LBS) with a subroutine which, 
for all cells, first determines the mass averaged velocity and thermal energy of all the particles 
in the cell. Then for each degree of freedom for each species of particle a new subroutine 
EQUIL (which executes the algorithm in Appendix A) is called in order to establish the 
new distribution of energy in the single degree of freedom amongst all the particles of the 
species. In so doing the mean velocity of each species is constrained to be exactly equal to 
the mass averaged velocity and the random energy stored in each degree of freedom is the 
same for each species and for each degree of freedom. Some implications of this are discussed 
in section (4). 

The input to the EQUIL subroutine is the number of particle velocities to be set N s , 
the required mean velocity of the single component of velocity u, and e' which is the total 
specific thermal energy stored in the single degree of freedom divided by the mass of one 
particle. Thus after the N s particle velocities have been set (for particles all having the same 
molecular mass) we have 

1 N 3 

« = ( 8 ) 

j = 1 

and 

1 N ‘ 

= o E ( u i “ ^) 2 - 

z j = i 

Note that u may be any one of the three components of molecular velocity (v x ,v y ,v z ) or 
a pseudo-velocity component representing a single degree of freedom for rotation, in which 
case u = 0. 

The required value of e' is determined before each call to EQUIL from 


( 9 ) 


€> — N $€■ d.o.f . j 


(10) 


where m s is the particle mass and e'd.o.f is the mean thermal energy per single degree of 
freedom given by 

e'd.o.f. = E'/v. (11) 
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In (11), E' is the total thermal energy of all particles in the cell and v is the total number 
of degrees of freedom. The latter is given by 

— IZO + M. (12) 

k= 1 

where £* is the number of fully excited (classical) degrees of freedom of the molecular struc- 
ture; thus £ = 0 for monatomic particles such as argon and helium and ( = 2 for diatomic 
species such as nitrogen and oxygen at room temperatures. 

In order to determine the total thermal energy E', the total translational energy E tr of the 
N particles in the cell must be divided into the bulk translational energy and the random 
thermal energy relative to the bulk motion. The thermal velocity* of the k-th particle is 
defined by 

v'* = v* - v est . (13) 

It follows from (3) and (13) that 

N N N 

Y m kK k = Y m k v 'y k = Y = °- ( 14 ) 

k= 1 fc=l k= 1 

The total translational energy of the N particles is given by 

E tr = \ Y m ^ 2 k = \ Y m k (< + V l k + V l k ) ■ ( 15 ) 

Z k=l Z k = 1 

By expressing the velocity of each particle in terms of the mass average velocity and thermal 
velocity (15) may be written as 

E tr = \ Y m k ((^ + V' xk ) 2 + (u, + V' yk ) 2 + (v 2 + <J 2 ) 

= \ Y ™k {vl + v] + v 2 z ) +\Y m « (v'l k + v'l k + v' 2 Zk ) 

* k= 1 C k= 1 

N N N 

+ V X Y +VyY + t», Y m k V Z k 

k=l k= 1 A:=l 

= - Mt^l^ + E' tr (16) 

where use has been made of the (14) and where 

£; = ix>*v'* 2 (i7) 

Z k = 1 

is the total translational thermal energy in the cell. The total energy of the particles is 
E t ot = E tr + E' rot , that is, it includes thermal energy stored in the rotation of the particles. 
The total thermal energy ( translation + rotation ) is given by 

E = E'„ + E' tM = Eut - (18) 

t Strictly speaking, v' is the estimated thermal velocity of a particle; that is, it is the velocity relative to 
the simulation estimate of the true gas velocity at that cell’s location. 
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4 Some comments on EPSM and DSMC 


It might be thought that another method of setting the equilibrium state within each cell 
could be used. Given that the estimates of the mean velocity and the temperature in each 
cell are known from (3) and (5), a set of equilibrium particle energies could be selected 
directly from the distribution (1) with v = v es< and T = T es t , in exactly the same way that 
an initial equilibrium state of a gas is set in standard DSMC. However in that case, the mean 
velocity of the finite sample so selected will not equal v est and the mean thermal energy of 
the finite sample of velocities will not correspond to T est ] in other words, momentum and 
energy will not be exactly conserved by the redistribution. Nanbu’s simulation method [8] 
is another direct simulation method that does not conserve momentum and energy exactly 
in collisions, and this seems to be the cause of the greater statistical scatter in the results 
compared with DSMC [9] [10]. It seems highly likely that the above procedure for setting an 
equilibrium state without preserving momentum and energy exactly would be unsatisfactory 
as well. 

The EPSM procedure does more than conserve energy as the new state in a cell is 
established. By setting each degree of freedom separately as described in section (3) the 
final state is slightly different from the state that would be established in DSMC after a 
large number of collisions had been calculated. In the DSMC equilibrium state the mean 
velocity of each species would be approximately, but not exactly, the same as the mass 
averaged velocity and the different degrees of freedom would not contain exactly the same 
thermal energy. There is some indication in the results that follow that EPSM shows less 
variation than DSMC between different runs (using a different sequence of random numbers) 
for the same flow. This may be related to the slightly less statistical variation that is allowed 
in EPSM when the velocity distribution function is represented by a finite sample of particle 
velocities. 

The EPSM algorithm could be used to set the initial state in DSMC in which case the 
initial state would contain an exact specified amount of energy and each degree of freedom 
would contain exactly the same total amount of thermal energy. It is possible that this might 
reduce the statistical scatter in DSMC calculations somewhat, particularly when data are 
obtained by averaging over an ensemble of different runs; in standard DSMC the nominally 
identical initial states contain slightly different total amounts of energy and momentum. It 
might also be possible, and useful, for something similar to be done to establish incoming 
fluxes at flow boundaries. In DSMC a selection of particles from an equilibrium gas flows 
into the simulation and this selection contains a statistically varying amount of energy and 
momentum at each step, even though the boundary condition is fixed; it might be preferable 
if exactly the same energy and momentum was carried in at each step (while the individual 
particle velocities varied between steps). 

5 Computational Efficiency of EPSM vs. DSMC 

The expected use of EPSM as part of a hybrid DSMC/EPSM code requires that EPSM be 
more efficient than DSMC for establishing equilibrium. In this section the computational 
efficiency of EPSM is investigated. In a hybrid code EPSM will only be used where the 
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standard DSMC procedures call for a large number of collisions to be calculated in a single 
time step. This can only happen if the normal requirement that At be small with respect 
to the local collision time is relaxed, at least in some part of the flow. If this criteria is not 
relaxed it follows that the average number of collisions per particle will be less than 1 and 
it would be inappropriate to use EPSM in such a situation unless it were known in advance 
that an equilibrium state existed in the cell. There are situations where large cell sizes are 
used in regions of flow where the flow gradients are small, and in that case the time step 
could be large relative to the local collision time. This is particularly likely to happen for 
multi-grid applications of DSMC. 

In order to compare computational efficiency, DSMC and EPSM have been used to 
simulate a gas at rest throughout the region between two stationary cylinders. The gas 
conditions were the same as the initial state considered in section (6) for Taylor- Couette 
flow. The expected result is that the flow remain in its initial state (in a statistical sense) 
for all values of the time step, and this was realized. The timing results are shown in figure 
(1). For values of At\f2RT\/ \\ = 10, for which the number of collisions per particle per 
time step is approximately 10, EPSM requires approximately 30% of the CPU time required 
by DSMC. Thus EPSM is an excellent alternative to DSMC in those cells in which many 
collisions must be calculated in one time step. It should also be noted that not a great deal 
of attention has yet been paid to considerations of computational efficiency in EPSM; there 
is probably more room for improvement. 


6 EPSM applied to Taylor- Couette Flow 

In order to demonstrate the use of EPSM an interesting vortical flow (Taylor- Couette flow) 
has been considered, even though the particular conditions considered are not necessarily 
those for which for which EPSM will be as accurate as DSMC. The classical case of Taylor- 
Couette flow occurs when a gas is contained between two coaxial cylinders of infinite length 
and the inner cylinder rotates while the outer cylinder is at rest. At low rotation speeds, a 
laminar flow is induced between the cylinders; this is an axisymmetric version of the familiar 
laminar Couette flow between two plates which move relative to each other. For higher 
rotational speeds a number of toroidal vortices will develop. At even higher rotational speeds 
the vortices can break down and the flow can become fully turbulent and three dimensional. 
Here axisymmetric flow only is considered; three velocity components are carried with the 
particles in the simulation (as is usual for DSMC) and the boundary condition at the surface 
of the rotating cylinder imparts a tangential velocity to the particles. Figure (2) shows the 
computational domain and an example of a time averaged EPSM flow with typical vortices. 

Stefanov and Cercignani [6] made a through study of this flow for a moderately rarefied 
gas and Bird [2] has repeated one of their calculations. Here both EPSM and DSMC sim- 
ulations have been performed. For DSMC the molecular model is a hard sphere with the 
molecular size appropriate to Argon at a temperature of 273 K [2]. The ordinary gas con- 
stant is R = kb/m = 208.2 J kg~ l K~ l . Argon has no rotational or vibrational energy and 
thus the ratio of specific heats is 7 = 5/3. 

The inner cylinder radius was 0.025 m and the outer cylinder radius was 0.05 m so that 
the gap between the cylinders was AR = 0.025 m. The simulations began with a uniform 
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gas at rest between the cylinders and the rotation of the inner cylinder began impulsively. 
The initial gas state was given by Xi = 273 K , and two values of initial density pi were used: 

1. pi = riim = 1.716 x 10“ 4 kg m~ 3 , P\ = 0.130 x 10 5 Pa , A* = 5 x 10 _4 m and 

2. p\ = nim = 2.574 x 10 -4 kg m~ 3 , Pi = 0.195 x 10 5 Pa , \i = 3.33 x 10 -4 ro. 

The tangential speed of the inner cylinder wall was 1011.33 m-s -1 . The inner and outer 
cylinder surface temperatures were constant at T w = Ti and diffuse reflection was assumed. 
In the axial direction the simulation region had a length of L — 4A R and was bounded by 
specularly reflecting surfaces at each end. 

For the given geometry the flow can be characterized by three non-dimensional parameters 

1. the Knudsen number based on the gap between the cylinders, Kn = Ai/A R = 0.02 or 
0.0133, 

2. the speed ratio S w = V w /\/2RJ\ = 3, where V w is the tangential speed of the inner 
cylinder and 

3. the wall to gas temperature ratio T w /Ti = 1. 

Different grids were used for the different Knudsen number flows; in both cases the cells 
were square with sides equal to Ai and the decoupling interval (time step) was set as At = 
~X/y/2RTi. Variable weighting factors were not used so the number of simulator particles 
remained constant throughout the simulation. It must be remembered that the EPSM 
simulations are independent of the collision rate (in effect the collision rate is infinite) and 
hence do not depend on the Knudsen number. The EPSM results may depend on the number 
of simulator particles used, the cell size and the decoupling interval. Some of these issues are 
addressed in section (7) which deals with the nature of the dissipation inherent in EPSM. 

The flow can be visualized with the aid of stream traces constructed from the velocity 
field in the x — r plane. Each drawing of these vortices used the same number and starting 
locations of the stream traces. Each figure shows the normalized time t = tV w / (2i rRi nner ) 
at the end of the sampling period. Except where explicitly stated in the caption the figures 
show results of averaging the flow over a single revolution of the inner cylinder. 

Comparison of DSMC and EPSM results gives an indirect method of determining how 
closely the DSMC calculations are approaching the equilibrium limit or how well the DSMC 
simulations represent the small departures from equilibrium which might be expected for 
small but non-zero Knudsen numbers. If the DSMC results were no different, or little 
different, from the EPSM results it would mean that the state in each cell was, in effect, in 
local equilibrium; it would then have to be considered whether this could be a valid result 
for this flow. If the results for EPSM and DSMC are different then it must be concluded 
that EPSM is not suitable for the case considered; this sort of information will be useful to 
determine the appropriate regime of application of EPSM in any hybrid DSMC/EPSM code. 
The question of the accuracy of the DSMC results depends on the usual requirements that 
the cell size and decoupling time step be small relative to the mean free path and collision 
time, respectively; these requirements are approximately, but not unambiguously, satisfied 
in these simulations. 
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It was found that EPSM required about 10% more CPU time than DSMC in these simu- 
lations. This is what can be expected according to figure (1) for the value of At used. This 
emphasizes the point that where the usual accuracy requirements can be met (or approxi- 
mately met), EPSM is, in theory, not only inferior (to an unknown extent) to DSMC but is 
also inefficient (at present) compared with DSMC; it is only where the usual requirements 
are relaxed, and the number of collisions per particle per cell becomes greater than 1 that 
EPSM is more efficient than DSMC and becomes strictly justified on the grounds that DSMC 
approaches the EPSM limit. 


6.1 Results for 200 x 50 grid. Kn = 0.02 

This grid is the same as that used in [2]. The flow was sampled at intervals of 3 At so an 
average particle in the undisturbed gas would move approximately 2 cell widths between 
samplings. In the time required for a single revolution of the inner cylinder the flow was 
sampled 52 times. 

The unsteady development of the flow calculated using EPSM is shown in figures (3), 
(4) and (5). At first there are four vortices which reduce to three after 10 revolutions and 
the sizes of the vortices vary with time. The development of the flow using DSMC is shown 
in figures (6), (7) and (8). DSMC gives a very different result from EPSM; there are five or 
six vortices rather than three. This is also different from the DSMC result of Bird [2], which 
contained three vortices, and of Stefanov and Cercignani [6], which contains four vortices. 
A DSMC flow is not exactly reproducible between runs and the same applies to an EPSM 
flow; both methods depend on the sequence of random numbers used in the simulation. 
Apparently DSMC is particularly sensitive to this effect [11] for this vortical flow. 

These calculations were repeated for both EPSM and DSMC using a different sequence of 
random numbers and the results after twenty revolutions are shown in figures (9) and (10); 
these can be compared with those shown in figures (5) and (8). Three vortices developed 
for all EPSM calculations made for this case, some of which continued for sixty revolutions. 
On the other hand, only four vortices developed for DSMC in this new run. The DSMC 
calculations used four sub-cells for each of the 50x200 cells in the grid§ and this gives an 
average of only three particles per sub-cell which seems a very small sample from which 
to be choosing collision partners; this small number of particles in each sub-cell may well 
be the source of the larger variation between runs for DSMC than for EPSM. The EPSM 
procedure was implemented over the entire cell; a better comparison between EPSM and 
DSMC could be made if EPSM were also implemented at the sub-cell level, which has not 
yet been done. It is worth noting that only about four particles per cell where used in 
Stefanov and Cercignani’s simulations [6] which also differ from the DSMC results here. 

6.1.1 Larger number of simulator particles 

This case ( Kn — 0.02, 50x200 grid) was re-calculated using twice the number of simulator 
particles for both EPSM and DSMC. Two different EPSM results are shown in figures (11) 
and (12). There is some difference between these results and each differs from the EPSM 

5 Collision partners are selected from the same sub-cell while the collision rate for the entire cell, not the 
individual sub-cells, is correctly simulated for the conditions averaged over the entire cell [2]. 
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result found with a smaller number of simulator particles for which there were three rather 
than the five or six vortices shown here. For DSMC there is a much bigger effect when the 
number of simulator particles is doubled; no vortices at all are formed (see figures (13) and 

(14))- 

On the assumptions 

1. that DSMC gives a better realization of this flow than EPSM does and 

2. that increasing the number of simulator particles increases the accuracy of the simula- 
tion 


(two assumptions that can hardly be doubted) there is a trend that the vortices become 
smaller (or disappear) as the accuracy increases. It might be that the simulation becomes less 
dissipative as it becomes more accurate and that more vortices develop for a less dissipative 
flow (in effect a higher Reynolds number flow) and eventually the vortices become unstable 
at sufficiently high Reynolds number. However, to exclude what seems a remote possibility, 
that the results for DSMC in this case correspond to a laminar flow with statistical noise, a 
detailed statistical analysis along the lines of that done in [6] would be necessary. 

The question of dissipation in EPSM is addressed in section (7); it is not immediately 
clear exactly how the number of simulator particles N to t is related to the dissipation in these 
methods, (although for partcile simulation methods the random noise 1 IsfNZt 12}). 


6.2 Results for 300 x 75 grid. Kn = 0.0133 

Both EPSM and DSMC simulations were repeated for a higher density flow and with a 
smaller cell size. The decoupling interval was such that At\/2R1 \/ remained at 0.67. The 
total number of particles was twice that for the 50x200 grid and the average number of 
particles per cell was 10.6, slightly less than the minimum on the 50x200 grid (11.9 particles 
per cell). The DSMC results would be expected to change for this different physical situation 
and the EPSM results would also be expected to change since, if the argument in section 
(7) is correct, the effective viscosity (effective mean free path) in EPSM depends on the cell 
size and decoupling interval. 

The unsteady development of the flow for EPSM is shown in figures (15), (16) and (17). 
This can be compared with the DSMC flow shown in figures (18), (19) and (20). For EPSM 
there is not much difference at t = 20 between the vortices in figure (17) for Kn = 0.0133 
and those in figure (9) for Kn = 0.02. and for DSMC there is not much difference between 
figure (20) for Kn = 0.0133 and figure (10) for Kn = 0.02. Of course a greater difference 
can be found by comparing with the low density runs of DSMC for which five vortices (or no 
vortices at all!) formed. If the purpose of this work were to investigate this Taylor- Couette 
flow in detail or to investigate fully the sensitivity of DSMC in this application it would 
obviously be desirable to repeat this low Knudsen DSMC calculation with a larger number 
of simulator particles; it could well be that the vortices disappear as they did for the higher 
Knudsen number case. 

While it is difficult to draw any firm conclusions about the Taylor- Couette flow from these 
results it is interesting to note that the DSMC results are different enough from the EPSM 
results (five vortices rather than three) to suggest that the velocity distribution function / 
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may be significantly different from equilibrium in these simulations even though the Knudsen 
number at 0.0133 is ‘small’. However it must be remembered that the product S w Iin is more 
relevant than Kn alone. The former may be expressed as 


S w I(n 


(2RTi) 


5 ' V,„ 


V W A 


(2 RTr) 


r/A-R 


(19) 


and is thus a ratio of 

1. a characteristic time between collisions t co u = A/ (2 RTi)* and a characteristic time of 
the macroscopic flow A R/V w or 


2. a characteristic distance V w t co u traveled (in the laboratory reference frame) by the 
molecules between collisions and a characteristic dimension of the macroscopic flow 
AFfl. 


While the usual definition of mean free path, the distance traveled between collisions mea- 
sured in a frame of reference moving with the local fluid velocity, makes it a convenient state 
property that combines the gas density, temperature and molecular size, it is not necessar- 
ily the most appropriate microscopic length scale for high speed flows, as should become 
apparent after writing the Boltzmann equation in non-dimensional form (see for example 
[12]). It seems unfortunate that the term ‘mean free path’ can not be used to refer to the 
distance between collisions measured in the whatever reference frame is relevant to the flow; 
if it could, the term ‘mean free path’ might be better applied to 

A = Vchartcoll (20) 

where i 

Acftor- — HiaX JV/Zutdi {‘2‘RTcharY 

is a characteristic speed based on either the macroscopic flow speed or a characteristic thermal 
speed of the molecules. Equation (20) may be written as 

A = max [1, S w ] A 

where A is the mean free path as conventionally defined. The mean collision time t co u or its 
inverse could serve as a state property which carries information about the gas state only, 
just as well as A can. 


7 Effective mean free path and dissipation in EPSM 

EPSM is an infinite collision rate limit of DSMC with a finite sample of simulator particles. 
The ‘infinite particle’ limit of EPSM is known as the Equilibrium Flux Method (EFM)H [5] 

^It is also the ratio of a characteristic shear stress in the flow, piV w / AR, and a characteristic pressure, 
piRTi and it may be written as ss M 2 /Re, where M is the Mach number and Re is the Reynolds number. 

II Not an exact limit, since some gradient information in EPSM carried by the possible non-uniform spatial 
distribution of particles in a cell is absent from EFM. Pullin’s EFM has has been rediscovered in [13] and 
given a new acronym (KFVS). Here I follow the convention that the naming rights belong to the originator 
of the method. 
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and has been used as an Euler flow solver in its first order [14] [15] [16] [17] and second 
order [18] [19] [20] formulations and as the basis of a Navier-Stokes solver [21] [22] [23]. 
The dissipation inherent in EFM has been demonstrated for a special case [15] and can be 
evaluated more generally [24]. It is worth reviewing those demonstrations since they shed 
light on the nature of the dissipation which arises in EPSM and, in some circumstances, in 
DSMC. 

If the decoupling interval (or time step At) in EPSM is set to a value larger than the 
time of flight across a typical cell then a typical particle will move to a new cell carrying 
with it momentum and energy which are representative of the conditions in the cell from 
which it came. Next its momentum and energy is transmitted to its new cell as it undergoes, 
in effect, an infinite number of collisions with all the other particles in the cell. Clearly 
the distance traveled between collisions in EPSM (and in DSMC) cannot be less than the 
distance traveled by a typical particle in time At. This corresponds to a minimum mean free 
path (as conventionally defined in the reference frame moving with the mean flow velocity) 
of approximately ( 8RT/n ) ? At. 

On the other hand, consider the case where At is small so that the chance that particles 
will cross more than one cell in any one time step is negligible. Consider two adjacent cells in 
an EPSM simulation as shown in figure (21). Let a local set of axes (n,p, q) be attached to 
the interface and let n point from the cell denoted (+) to the cell denoted (— ). Let v n , u p , v q 
be the components of molecular velocity with respect to the local axes. All particles having 
Cn > 0 which cross the interface are selected from the equilibrium distribution established 
in the (+) cell; all particles crossing the interface having c n < 0 are selected from a different 
equilibrium distribution Jq. As a result the distribution of Cn of those particles which cross 
the interface is discontinuous at c„ = 0 as shown in the figure. 

It is to be expected that the fluxes derived from this (discontinuous, non-equilibrium) 
distribution of velocities of those particles crossing the interface will contain parts corre- 
sponding not only to the continuum Euler fluxes (for a hypothetical equilibrium state at the 
interface) but extra ‘dissipative’ fluxes. Note that figure (21) shows a hypothetical equilib- 
rium distribution derived from the non-equilibrium distribution; the full set of dissipative 
fluxes for EPSM (relative to this hypothetical equilibrium distribution) is given in [24] and 
these results could be used to derive an effective mean free path for EPSM. Here, to illustrate 
the point, a simple example of the exchange of fluxes between cells in EPSM is considered. 

7.1 Fluxes in EPSM 

Let Q be one of the molecular properties m,mv or m (|v 2 + e rot j (mass, momentum and 
energy) which are transported as the particles move. The net flux (/unit area/unit time) 
of Q across the interface can be split into two parts: a forward flux Fq which flows in the 
direction of n and a backward flux Fq which flows in the opposite direction. The total flux 
is 

Fq = F+ + Fq ( 21 ) 

where the two parts of the fluxes can be obtained from standard kinetic theory ( e.g . [2]) as 

/ oo yoo yoo 

/ / n + Qv n dv n dv p dv q (22) 

■oo J — oo JO 
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and 


/ oo roo rO 

I / ^ fo Qv n dv n dvpdvq. (23) 

-oo J — oo J — oo 

In (22) and (23), n = p/m is the number density of particles and m is the mass of one 
particle. 

Now consider the idealized situation in figure (22) which shows adjacent cells in a section 
of the flow where there the gradient dv p /dx n ^ 0 To evaluate the net transfer of p-momentum 
across the interface we need to put Q = mcp and perform the integrations in (22) and (23). 
With the form of and known from (1) the forward and backward fluxes of p-momentum 


are s 

FL, = W'Vttft’? ± £V (2 RT^vp 

where 

IV ± = i[l±erf(4)], 

(24) 

(25) 

B± - ± 2,| eXP [ W)]‘ 

(26) 

and 


4 = 4/ ■ 

(27) 

If we further assume there are no gradients other than dv p /dx n , we have p + = p = p, 

= v n = v+ = v q =0 and T + = T = T, and the net flux (/unit 
p-momentum across the interface is 

area/unit time) of 

F„„ =p(RT/2*)i („+ - v;) . 

(28) 

This net flux is the shear stress r acting at the interface and, since 
approximated by Ax n dv p /dx n where Ax n is the cell size, we have 

{y* -v~J may be 

r = p(RT/2ir)* Ax n dv p /dx n . 

(29) 


Comparing this to the stress-strain relationship for plane Couette flow, r = p \dv p /d\ x n , we 
get the effective viscosity for EPSM as 

p' = (RT/2tt)i pAx n . (30) 

Comparing this to the usual relation between mean free path and viscosity 

p ~ -p (8RT/^ A 

we obtain an effective mean free path in EPSM of 

A' » Ax n /2. (31) 
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Combining the results for At > Ax/v and At < Ax/v where Ax is a characteristic cell 
size and v is the magnitude of the local mean particle velocity (mean fluid velocity) we can 
say that there is an effective mean free path in EPSM given by 


A epsm ~ max 


^ Ax, (8RT '/tt) * At 


1 

— max 
£ 


I c 

7T2 O . 


Ax 


(32) 


where C = vAt/Ax is the CFL number and S = v/ ( 2RT ) 2 is the local speed ratio. We 
could also say that the effective viscosity in EPSM is given by 


ft epsm ~ max 


i ±£ 

*JSJ 


IRT/2-p P Ax. 


(33) 


7.2 Implications for DSMC 

The above analysis applies approximately to DSMC and becomes more exact as the distribu- 
tion function approaches an equilibrium form. Hence the effective mean free path in DSMC 
would be given by 

A DSMC = max [ 0 (Ax„) , A] 

where A is the local mean free path for the molecular collision model being used. This is 
merely another way of stating the theoretical requirement that the cell size in DSMC should 
be less than A. It is worth noting, however, that the effective viscosity in DSMC will be 
given by an expression similar to (30) and not by the theoretical viscosity of the molecular 
collision model, when the cell size is larger than A, as may be the case when it is based on a 
(larger) local characteristic length derived from flow gradients. In other words, the collision 
model is largely irrelevant when the cell size is large, and in that case EPSM is probably 
little different from DSMC. 

It is obvious that DSMC and EPSM should give similar results when the number of 
collisions per particle per time step in DSMC is large, which is likely to be the case when 
cell sizes, and hence the decoupling interval, are large. This does not necessarily imply an 
inadequacy in DSMC or EPSM, since when flow gradients are small and large cell sizes 
are used, the local flow state is most likely to be very close to equilibrium. An analogy in 
continuum fluid mechanics is that the Euler equations (rather than the full Navier-Stokes 
equations) can give reasonable results in regions of the flow where gradients are small. 


8 Conclusions 

The equilibrium particle simulation method (EPSM) [5] has been extended to handle flows 
with multiple species and multiple (classical) degrees of freedom of molecular structure. 
EPSM is now a suitable method to be used as the high density part of a hybrid particle 
simulation method for flows which contain regions of rarefaction as well as regions of near 
continuum flow. EPSM cannot be expected to be used as an alternative to DSMC for all 
flows for which DSMC might be used. It is restricted theoretically to those cases (which 
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may be only a small fraction of the cells in a hybrid DSMC/EPSM simulation) where the 
number of collisions per particle per time step is large and hence the distribution function 
in a cell can be expected to take an equilibrium form. In that case it has been found that 
EPSM is very efficient compared to DSMC. 

It was not the purpose of this work to develop a hybrid DSMC/EPSM code but merely 
to demonstrate the use of EPSM. To this end, the Taylor Couette flow of a slightly rarefied 
gas (An = 0.02 and 0.0133) at high speed ( S w = 3) has been calculated using EPSM 
and DSMC. Differences between DSMC and EPSM results were found and it follows that 
even at the low Knudsen number of 0.0133 the simulated flow state is non-equilibrium. 
What is more important than the Knudsen number alone for determining the degree of 
departure from equilibrium in high speed flows is the product S w Kn which is the ratio of 
the characteristic time between collisions and a characteristic time of the macroscopic flow; 
in the flows considered here this ratio has a value of almost 4% which is apparently enough 
to produce a significant departure from equilibrium in the simulation. It doesn’t appear safe, 
based on the variation between different DSMC results, to make strong conclusions about 
Taylor- Couette for the conditions considered here. Great computational effort would need 
to be expended (particularly on larger numbers of simulator particles but probably also on 
smaller cell sizes). 

EPSM was found to be less sensitive to a change in the number of simulator particles used 
than was DSMC. There was some indication that EPSM was more dissipative than DSMC, 
possibly because EPSM was implemented at the cell level, rather than the sub-cell level for 
which collisions are calculated in DSMC. The question of the accuracy of EPSM in the near 
continuum regime has been addressed theoretically by considering the nature and magnitude 
of the effective dissipation and mean free path in EPSM. This analysis serves to emphasize 
some possible problems in using particle simulation methods (both EPSM and DSMC) in 
near continuum flows where the usual requirement that the cell size and decoupling interval 
both be small is relaxed. It that case the dissipation in the simulation methods could depend 
more on the cell size than on the collision model being used; in this respect the simulation 
methods become more like continuum solution methods ( e.g . [15]) which contain an inherent 
‘artificial’ or ‘grid’ dissipation. 
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Appendices 


A The EPSM (equilibrium) algorithm 

The core of the EPSM code is the the EQUIL subroutine which establishes equilibrium values 
for one degree of freedom for a finite sample of particles. The derivation of the method is 
given in [7] and a description of the algorithm is given in [5]. A more detailed pseudo-code 
description is given below. Each call of TZ returns a new random fraction between 0 and 
1 selected from a uniform distribution. The input to the EQUIL subroutine is the number 
of velocity components to be set N, and the required values of u and the variance e! . The 
output is the set of values u^j = 1 , N which are selected from a normal distribution and 
which satisfy the constraints u = Yl u j/N and e' = Yl( u j — u) 2 /2 (see (8) and (9)). The 
algorithm may be easier to follow using from the following diagram which illustrates, for 
particular values of N, the order of generation and interdependence of the various terms in 
the algorithm. 


N = 10 _ N = 9 

e' u e' u 



Interdependence of the values ei,Vj and Uk for N = 10 and N = 9. 
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if 1V = 

2 then 

Hi 

<— u ± \fe' (equal probability for ± 

u 2 

<— 2 u — u\ 

else if = 3 then 

e 

<- 2 irR 

v 2 

<— \/2e'cos0 

V Z 

<— V^^sin# 

U X 

<— u — \/2 V 2 /\/3 

U 2 

«- U\ + (\/3h2 - W 3 ) / \/2 

Uz 

n 2 + v/2^3 

else if TV = 4 then 

Ti 

7e 2 

^2 

-e'(l-r x ) 

e 

«- 2tt^ 

v 2 

<— -v/2e 2 cos 6 

Vz 

<— \/2e 2 sin 6 

Ui 

<r- u — \/3v 2 /\/iv 2 

U 2 

<— ui + (\/4 u 2 — \/2vzj /\/3 

e 4 

^eTx 

v 4 

<— ±\/2e7 (equal probability for ± ) 

U Z 

*— u 2 4- (\/3 V 3 — v 4 j / \/2 

U 4 

< — U 3 + \/2 u 4 


else if TV > 4 then 


if N is even then 

k max *— N /* maximum index of */ 

for m = 1, . . . , k max / 2 — 1; step of 1 

1 

T m < — 'J^k-max I"*— m—0.5 

end for 

else if N is odd then 

kmax «— N — 1 /* maximum index of tk */ 

for m = 1, . . . , k max /2 — 1; step of 1 
T m R k max/i-m 
end for 
end if 

e 2 «— e'(l — Ti) 

9 «- 2ttR 

v 2 <— \/ 2 e 2 Cos^ 

^3 «— \Z5e7sin0 

«x <— u — \/TV — 1 v 2 /y/N 

u 2 *- ui + (s/N v 2 - \JN -2v^j I'JN - 1 

continued next page .... 
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continue general case for N > 4 .... 

P 1 

/* initialize the product P = flmli 1 */ 

for k = 4, . 

■ ■ , k max — 2; step of 2 /*last value of e^, k = k max , is done later 

P 

T 

'to 

Ck 

e ' (l ~ Tk/2 j P 

e 

<— 2nlZ 

Vk 

<— \/2e^ cos 6 

Vk+1 

<r— yjltk sin 9 

Uk-1 

<— Uk- 2 + f - \/-/V + 3 — k Vk-i — VN 4- 1 — k Vk) /y/N + 2 — k 

Uk 

<— u fc _i + [y/N + 2 - - v/iV - /V7V + 1 - fc 

end for 
P *- 

Tkmxx/2-lP 

e kmax * 

e'P 

if N is even then 

v N 

<— ± y /2ek Jnax (equal probability for ±) 

else if N is 

odd then 

e 

«- 27r7e 

Vn-1 

«- V 2e kmax COS 9 

VN 

*- V 2e kmax sin 9 

UN-2 

<— un — 3 + (y/i VN- 2 - >/2 u^-i) / ^(3) 


end if 

UN-1 Un - 2 + (v/3 VN-l ~ Vn) / \/2 

un <— un-i + \/2 vn 
end if /*end of general case N > 4 */ 
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B Some coding errors in DSMC2A 

A number of modifications have been made to the subroutine M0VE2A which calculates 
the movement of the particles in the axisymmetric flowfield. The most important change 
concerns the procedures to be followed after a collision with a reflecting surface has been 
detected. The particle is moved to the collision point by calling subroutine AIFR. The input 
to AIFR is the initial radial coordinate r (which corresponds to the initial y coordinate in 
3D space), the initial radial and tangential velocity components (which correspond to the v y 
and v z velocity components in 3D space), and the position coordinate changes Ay and A z 
which occur during the time of flight up until collision with the surface. The output from 
AIFR is the new radial coordinate and the new radial and tangential velocity components 
(there is a continuous alteration of the radial and tangential coordinate directions as the 
particle moves trough 3D space with constant velocity). 

The first apparent error in M0VE2A occurs in the call to AIFR when a collision is 
detected with a surface which is normal to the x axis; AIFR is called with increments Ay 
and Az which are appropriate to the time of flight after rather than before the collision. The 
second, and more serious error, concerns the checking of the weighting factors after surface 
collisions to determine whether the particle (numbered N) should be deleted after collision 
or duplicate particles should be created. In that the particle numbered N is moved for the 
time of flight remaining after collision on return from WEIGHT, the code implicitly assumes 
that the particle was not deleted. If in fact the particle was deleted inside WEIGHT the 
value of N will have been decremented by 1 and so a lower numbered particle, which will 
have been moved previously in the given time step, will be moved again. 

This, perhaps small, error can lead to a more serious effect. If the first numbered par- 
ticle N = 1 happens to collide with a reflecting surface, and happens to be deleted inside 
WEIGHT, N will have the value 0 after the return from WEIGHT. The code will attempt 
to move the non-existent particle N = 0; this leads to references to memory locations which 
are ‘out of bounds’ of various arrays. To avoid this some extra few lines of code have been 
added which detect whether the particle was deleted inside WEIGHT; if so the code moves 
on to the next particle** 

Another very small error in the conservation of mass was detected. It was noticed that 
for the closed flow volume between the two rotating cylinders of the Taylor- Couette flow 
the total number of particles in the simulation was not constant, even though no weighting 
factors were used; there was a tiny percentage of simulator particles which were removed 
from the simulation. To overcome this, a number of small changes were made to the surface 
collision detection logic. These changes were designed to catch some singular cases when 
the particle happened to pass exactly (to machine precision) through the end point of a 
reflecting surface and when the particle finished its flight exactly on the surface. A similar 
problem occured when a particle would finish its flight just beyond a surface; it one part 
of the code this would be detected as a valid collision (which it is) but in a subsequent 
part of the code (because of rounding errors) this case was rejected as a collision and the 

“The general axisymmetric code in use at Imperial College for many years [25] [26], which was derived 
from an early version of the codes now designated G2 in [2], checked weighting factors only when the final 
position of the particle was found. Implementing this strategy in the DSMC2A code is another possible way 
of avoiding the problem. 
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particle was left in its position just beyond the surface. Since the surface lay along the 
boundary of the simulation region the particle was subsequently deleted. A change to 8- 
byte (double precision) arithmetic when solving a quadratic equation which detects surface 
collisions seems to have eliminated the known occurrences of this elfect. 

The detection of multiple reflections from surfaces and boundaries in an axisymmetric 
flow is a complicated task which was addressed in a report [27] on the 2-D DSMC code 
which was in use at Imperial College, London. The DSMC2A code has not been checked 
against all the possible errors mentioned in [27] and there may still be some remaining 
errors (or some new ones may have been introduced). However, after extensive testing in an 
enclosed flow with no weighting factors, it has become reasonably clear that any remaining 
ways for particles to be erroneously deleted are fewer than before. The original frequency of 
occurrence was small and is smaller now. 

An apparent typographical error in the main program (BFND for FND) where a previ- 
ously computed flow state was read from a file, has been corrected. The extent to which the 
errors mentioned here may be repeated in other DSMC codes supplied with [2] has not been 
investigated. 
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Figure 1: CPU time (SPARCstation 10) required to move 119,124 simulator particles 156 
times ( tmax — 156 At ), sample the flowfield 52 times, and write the final flow state to disk. 
Gas at rest between two coaxial stationary cylinders. For DSMC the non-dimensional time 
step shown on the x-axis is approximately equal to the number of collisions per particle per 
time step. 
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Gas between cylinders 


Computational domain 


Inner cylinder rotates 



diffusely reflectina surface T... = T, 



diffusely reflecting surface T w = T, 


Figure 2: Computational domain and time averaged simulated flow showing toroidal vortices. 
EPSM. 200x50 grid. N to t = 119, 124. t = 41—60 revolutions. 


24 


symmetry b.c. 


0.00 0.02 0.04 0.06 0.08 0.10 

Figure 3: Kn = 0.02. 200x50 grid. EPSM. N to t = 119,124. t = 5 revolutions. 
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Figure 4: Kn = 0.02. 200x50 grid. EPSM. N to t = 119, 124. t = 10 revolutions. 
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Figure 9: Repeat run with different random numbers. Kn = 0.02. 200x50 grid. EPSM. 
N tot = 119, 124. t = 20 revolutions. 
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Figure 11: Average of 23 simulator particles/cell. Kn = 0.02. 200x50 grid. EPSM. N-, 
238, 248. t = 10 revolutions. 
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Figure 12: Average of 23 simulator particles/cell. Kn = 0.02. 200x50 grid. EPSM. N t 
238,248. t = 10 revolutions. Repeat run with different random numbers. 
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Figure 13: Average of 23 simulator particles/cell. Kn — 0.02. 200x50 grid. DSMC. N to t 
238,248. t = 10 revolutions. 
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Figure 14: Average of 23 simulator particles/cell. Kn = 0.02. 200x50 grid. DSMC. N tot 
238,248. t — 10 revolutions. Repeat run with different random numbers. 
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Figure 18: Kn = 0.0133. 300x75 grid. DSMC. N to t = 238,248. t = 5 revolutions. 
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Figure 19: Kn = 0.0133. 300x75 grid. DSMC. N to t = 238,248. £ = 10 revolutions. 
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Figure 20: Kn = 0.0133. 300x75 grid. DSMC. N to t = 238,248. i = 20 revolutions. 
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Molecular velocity normal to cell interface v n /(2RT + ) 1/2 

Figure 21: Distribution of Cn for particles crossing the interface between two cells. Solid 
lines: distribution in EPSM; dashed line: hypothetical equilibrium distribution derived from 
the EPSM distribution 
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Figure 22: Two cells exchange momentum giving rise to an effective shear stress between 
them. 
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